In [1]:
import numpy as np
from numpy import random
import math
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline 

#matrix factorization packages
from sklearn.decomposition import FastICA, NMF
from pysptools import eea

In [2]:
#draw 559 samples from a normal Gaussian distribution 
random_noise=np.random.normal(loc=1, size=(559))

Define parameters for normal distribution


In [3]:
mu=400; variance=250; sigma=math.sqrt(variance); x=np.linspace(0,559,559)
#mu2=0; variance2=1; sigma2=math.sqrt(variance2); x2=np.linspace(mu2-5*sigma2, mu2+5*sigma2, 100)
mu2=250; variance2=100; sigma2=math.sqrt(variance2); x2=np.linspace(0, 559, 559)

Create normal distributions. These will be constituent spectra


In [4]:
#segmentA=mlab.normpdf(x, mu, sigma).reshape(559,1)*.2
#segmentB=mlab.normpdf(x2, mu2, sigma2).reshape(100,1)*.5
constituentA=np.zeros((559,1))
#constituentA[20:80]=segmentA
constituentA=mlab.normpdf(x, mu, sigma).reshape(559,1)*.5+mlab.normpdf(x, mu-3*sigma, sigma).reshape(559,1)*.25
constituentB=np.zeros((559,1))
#constituentB[400:500]=segmentB
constituentB=mlab.normpdf(x2, mu2, sigma2).reshape(559,1)*.5

In [5]:
plt.plot(constituentA)
plt.plot(constituentB)


Out[5]:
[<matplotlib.lines.Line2D at 0x4912630>]

In [6]:
test=abs((random_noise.reshape(559,1)/1000)+constituentB)
plt.plot(test)


Out[6]:
[<matplotlib.lines.Line2D at 0x194a8ba8>]

Make hyperspectral matrix


In [7]:
hyper=np.zeros((64,64,559))
for x in range(hyper.shape[0]):
    for y in range(hyper.shape[1]):
        if x%2==0:
            hyper[x,y,:]=constituentA.reshape(559)+abs(random_noise.reshape(559)/1000)
        else:
            hyper[x,y]=hyper[x,y]
        if y%3==0:
            hyper[x,y]+=constituentB.reshape(559)+abs(random_noise.reshape(559)/1000)
        else:
            hyper[x,y]=hyper[x,y]

#totally random spectra -- independent components
hyper2=np.zeros((64,64,559))
for x in range(hyper.shape[0]):
    for y in range(hyper.shape[1]):
            hyper2[x,y,:]=random.rand(1)*constituentA.reshape(559)+abs(random_noise.reshape(559)/1000)\
            +random.rand(1)*constituentB.reshape(559)+abs(random_noise.reshape(559)/1000)

#random, dependent components
hyper3=np.zeros((64,64,559))
for x in range(hyper.shape[0]):
    for y in range(hyper.shape[1]):
        frac=random.rand(1)
        hyper3[x,y,:]=frac*constituentA.reshape(559)+abs(random_noise.reshape(559)/1000)\
        +(1-frac)*constituentB.reshape(559)+abs(random_noise.reshape(559)/1000)
        
#pure regions with random spectra surrounding
hyper4=np.zeros((64,64,559))
hyper4[12:20, 45:50]=constituentA.reshape(559)
hyper4[42:62, 2:22]=constituentB.reshape(559)
for x in range(hyper.shape[0]):
    for y in range(hyper.shape[1]):
        if hyper4[x,y,:].all()==0:
             hyper4[x,y,:]=random.rand(1)*constituentA.reshape(559)+abs(random_noise.reshape(559)/1000)\
            +random.rand(1)*constituentB.reshape(559)+abs(random_noise.reshape(559)/1000)
            
#pure regions with dependent spectra surrounding
hyper5=np.zeros((64,64,559))
hyper5[12:20, 45:50]=constituentA.reshape(559)
hyper5[42:62, 2:22]=constituentB.reshape(559)
for x in range(hyper.shape[0]):
    for y in range(hyper.shape[1]):
        if hyper5[x,y,:].all()==0:
            frac=random.rand(1)
            hyper5[x,y:]=frac*constituentA.reshape(559)+abs(random_noise.reshape(559)/1000)\
            +(1-frac)*constituentB.reshape(559)+abs(random_noise.reshape(559)/1000)

In [8]:
plt.plot(constituentA)


Out[8]:
[<matplotlib.lines.Line2D at 0x1999d1d0>]

In [9]:
hyper_list=hyper4.reshape(4096,559)

ICA


In [10]:
ica = FastICA(n_components=2, max_iter=400)
#reconstruct signals
S = ica.fit_transform(hyper_list)
#mixing matrix
A = ica.mixing_ 
ICA_comp=ica.components_

In [11]:
plt.plot(ICA_comp[0,:])
plt.plot(ICA_comp[1,:])


Out[11]:
[<matplotlib.lines.Line2D at 0x1976ae10>]

Nonnegative Matrix Factorization


In [45]:
NMF_model=NMF(n_components=2)
NMF_fit=NMF_model.fit_transform(hyper_list)
NMF_comp=NMF_model.components_

In [46]:
plt.plot(NMF_comp[0,:])
plt.plot(NMF_comp[1,:])


Out[46]:
[<matplotlib.lines.Line2D at 0x2526a400>]

Automatic Target Generation


In [ ]:
ATG=eea.ATGP()
ATG_comp=ATG.extract(hyper4, q=2)

In [ ]:
plt.plot(ATG_comp[0,:])

In [ ]:
plt.plot(ATG_comp[1,:])

Pixel Purity Index


In [125]:
PPI=eea.PPI()
PPI_comp=PPI.extract(hyper,q=2)

In [126]:
plt.plot(PPI_comp[0,:])


Out[126]:
[<matplotlib.lines.Line2D at 0x2ad9fcc0>]

In [127]:
plt.plot(PPI_comp[1,:])


Out[127]:
[<matplotlib.lines.Line2D at 0x2af80668>]

Fast Iterative Pixel Purity Index (FIPPI)


In [284]:
FIPPI=eea.FIPPI()
FIPPI_comp=FIPPI.extract(hyper4, q=2)

In [285]:
plt.plot(FIPPI_comp[0,:])


Out[285]:
[<matplotlib.lines.Line2D at 0x39b976a0>]

In [286]:
plt.plot(FIPPI_comp[1,:])


Out[286]:
[<matplotlib.lines.Line2D at 0x39e73be0>]

N-FINDR


In [138]:
NFINDR=eea.NFINDR()
NFINDR_comp=NFINDR.extract(hyper, q=5)

In [139]:
plt.plot(NFINDR_comp[0,:])


Out[139]:
[<matplotlib.lines.Line2D at 0x2bccddd8>]

In [140]:
plt.plot(NFINDR_comp[1,:])


Out[140]:
[<matplotlib.lines.Line2D at 0x2c0b2d30>]

In [ ]: